18 Production model - river 4 (IS)

Simple CJS model integrated with a growth in weight model to get phi, p, and growth estimates to develop a production model.

In all the models below, 1 = not observed and 2 = observed in the input encounter data.
Encounter data are available here in the eh.csv file. Weight data are in weight.csv

Using model #3 from modelsCMR_Growth_NN_OB.qmd as a staring point for the models, but adapting the model by
Prob not 1) Extending AgeInSamples from 1-11 to 1-x to allow bigger fish to be present for the production estimates. Prob not 2) Loop over first[i]:lastAIS so fish have the chance to survive to large size.
3) Add sample random effect structured by cohort to allow cohort effects on growth.

18.1 By Cohort, modelNum 3: phi(i,t) * g(i,t), p(i,t) model

Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)

Probability of survival (phi) model structure:

      logit(phi[i,t]) <- 
        phiInt[i,t] + 
        phiBeta[1,i,t] * weight[i,t] + 
        phiBeta[2,i,t] * weight[i,t] * weight[i,t] +
        phiBetaCohort[cohort[i]]

Probability of capture (p) model structure:

p(t,i) <- pInt(t,i)

Growth rate (gr) model structure:

gr(t,i) <- grInt(t) + cohort(i)

Survival/growth rate interaction:

Multiplicative

Model code is in ./models/production/modelProduction_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/modelProduction_OB_makeFile.R
Functions for this qmd file in ./models/qmdProduction_OB_functionsToSource.R

Model runs:

18.1.1 How many ageInSamples to include?

Code

all <- tar_read(cdWB_electro_target)
table(all %>% filter(river == "wb obear")|> dplyr::select(ageInSamples))

ageInSamples
   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
1308 1890 1007  837 1486  875  446  298  444  159  155   57   47   17   16    3 
  16   18 
   3    3

Code

cohorts <- 2002:2014

18.1.2 Retrieve model results

Code

library(targets)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
# 
# 

# To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
source('./models/production/modelProduction_OB_functionsToSource.R')
source('./models/production/qmdProduction_OB_functionsToSource.R')

modelNum <- 3 # phi * growth

18.1.3 Model code

Code

# all cohorts have the same code. show here
load('./models/production/runsOut/production_OB_model_3_2002__mostRecent.RData')
d$modelCode

{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]], 
            sd = 2)
        weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = weightSD)
        for (t in (first[i]):(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]/90
            weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1], 
                sd = weightSD)
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
        }
    }
    for (t in 1:(T - 1)) {
        grIntT[t] ~ dnorm(0, sd = 1000)
    }
    delta[1] <- 1
    delta[2] <- 0
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gamma[1, 1, t, i] <- phi[i, t]
            gamma[1, 2, t, i] <- 1 - phi[i, t]
            gamma[2, 1, t, i] <- 0
            gamma[2, 2, t, i] <- 1
        }
    }
    for (t in 1:(T - 1)) {
        p[t] ~ dunif(0, 1)
        omega[1, 1, t] <- 1 - p[t]
        omega[1, 2, t] <- p[t]
        omega[2, 1, t] <- 1
        omega[2, 2, t] <- 0
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i, 
                t] * weight[i, t] + phiBeta[2, i, t] * weight[i, 
                t] * weight[i, t]
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
            for (b in 1:2) {
                phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
            }
        }
    }
    for (t in 1:(T - 1)) {
        phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
        phiIntTOut[t] <- ilogit(phiIntT[t])
        phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
        phiBetaT[2, t] ~ dnorm(0, sd = 1/0.667)
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            weightZ01[i, t] <- weight[i, t] * z[i, t]
        }
    }
    for (i in 1:N) {
        z[i, first[i]] ~ dcat(delta[1:2])
        for (j in first[i]:(last[i] - 1)) {
            z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
            yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
        }
    }
}

Code

getSummaryByCohort <- function(cohort, modelNum = 3) {
  # Load the data
  load(paste0('./models/production/runsOut/production_OB_model', '_', modelNum, '_', cohort, '__mostRecent.RData')) 
  
  # Input data
  eh <- tar_read_raw(paste0('eh_OB_', cohort, '_target'))
  
  get_MCMCTrace_OB(d)
  
  MCMCplot(object = d$mcmc, params = c("phiIntTOut"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("phiBetaT"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("p"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("grIntT"), ref_ovl = TRUE)
  
  # Create the summary list
  summary_list <- list(
    eh = eh,
    stats = kable(as_tibble(d$runData), caption = "Run statistics"),
    runTime = paste0('Run time = ', round(d$runTime, 3), ' ', attr(d$runTime, "units")),
    #traces = get_MCMCTrace_OB(d),
    summaryMod3_tt_growth = MCMCsummary(object = d$mcmc, params = c("phiIntTOut", "p", "grIntT", "phiBetaT"), round = 3) %>%
      rownames_to_column(var = "var"),
    summary_OB = get_summary_OB(d, eh)
  )
  
  return(summary_list)
}

#s=summary_OB[, order(summary_OB$ind, summary_OB$t)]

18.1.4 Summaries by cohort

Code

s2002 = getSummaryByCohort(2002, modelNum)

Code

s2003 = getSummaryByCohort(2003, modelNum)

Code

s2004 = getSummaryByCohort(2004, modelNum)

Code

s2005 = getSummaryByCohort(2005, modelNum)

Code

s2006 = getSummaryByCohort(2006, modelNum)

Code

s2007 = getSummaryByCohort(2007, modelNum)

Code

s2008 = getSummaryByCohort(2008, modelNum)

Code

s2009 = getSummaryByCohort(2009, modelNum)

Code

s2010 = getSummaryByCohort(2010, modelNum)

Code

s2011 = getSummaryByCohort(2011, modelNum)

Code

s2012 = getSummaryByCohort(2012, modelNum)

Code

s2013 = getSummaryByCohort(2013, modelNum)

Code

s2014 = getSummaryByCohort(2014, modelNum)

18.1.5 Combine cohort summaries

Code

sAll <- tibble(.rows = 0) |>
  add_column(!!!s2002$summary_OB[0,])

# Fill tibble
for (cohort in cohorts) {
  summary_obj <- get(paste0("s", cohort))$summary_OB
  summary_obj$cohort <- cohort
  sAll <- bind_rows(sAll, summary_obj)
}

sAll$cohort_ind <- paste0(sAll$cohort, "_", sAll$ind)

Code

#ojs_define(numTags_OJS_mod3 = dim(s2002$eh$tags)[1]) # all cohorts have the same eh
ojs_define(summary_OB_OJS_all = sAll)
#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))

18.1.6 Select cohort(s):

Code

summary_OB_OJS_all_T = transpose(summary_OB_OJS_all)

Code

import { range } from "d3";

viewof selected_cohort = Inputs.select(d3.range(2002, 2015), {
  label: "Which cohort?",
  value: 2002,
  step: 1,
  multiple: true
})

Code

summary_OB_OJS_all_T_selected_cohort = summary_OB_OJS_all_T.filter((d) =>
  (selected_cohort.includes(d.cohort))
)

18.1.7 Select one or more individuals:

Code

uniqueInds = [...new Set(summary_OB_OJS_all_T_selected_cohort.map(d => d.cohort_ind))].sort();

viewof selectInd_mod3 = Inputs.select(uniqueInds, {
  label: "Which fish?",
  value: 1,
  step: 1,
  multiple: true
})


summary_OB_OJS_all_T_selected_cohort_ind = summary_OB_OJS_all_T_selected_cohort.filter((d) =>
  (selectInd_mod3.includes(d.cohort_ind))
)

Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

18.1.8 Plot survival

Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code

Plot.plot({
    width: width,
    height: 350,
    inset: 10,
    color: {
      scheme: "greys"
    },
    x: { label: "Age/season combination" },
    y: { label: "Probability of survival" },
    marks: [
      Plot.frame(),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind, {
        x: "t",
        y: "pSurv"
      }),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind, {
        x: "t",
        y: "enc8",
        fill: "orange"
      }),
      Plot.line(summary_OB_OJS_all_T_selected_cohort_ind, {
        x: "t",
        y: "pSurv"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind, {
        x: "first",
        y: 1,
        "stroke": "green"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind, {
        x: "last",
        y: 1,
        "stroke": "red"
      })
      ,
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind, 
                Plot.selectFirst({
                  x: 11,
                  y: 1,
                  text: d => d.cohort
                })
      )
    ],
    facet: {
      data: summary_OB_OJS_all_T_selected_cohort_ind,
      x: "ind"
    }
  })

Code

Plot.plot({
    width: width,
    height: 350,
    inset: 10,
    color: {
      scheme: "greys"
    },
    x: { label: "Age/season combination" },
    y: { label: "Standardized body mass (mg)" },
    marks: [
      Plot.frame(),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind, {
        x: "t",
        y: "mean_weight"
      }),
      Plot.line(summary_OB_OJS_all_T_selected_cohort_ind, {
        x: "t",
        y: "mean_weight"
      }),
      <!-- Plot.dot(d, { -->
          <!--   x: "t", -->
          <!--   y: "mean_gr", -->
          <!--   fill: "green" -->
          <!-- }), -->
        Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind, {
          x: "t",
          y: "weightDATA_std",
          fill: "orange"
        }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind, {
        x: "first",
        y: 3,
        "stroke": "green"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind, {
        x: "last",
        y: 3,
        "stroke": "red"
      }),
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind,
                Plot.selectFirst({
                  x: 1,
                  y: 4,
                  frameAnchor: "top-left",
                  text: (d) => "Cohort = " + d.cohort
                })
      ),
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind,
                Plot.selectFirst({
                  x: 1,
                  y: 3.5,
                  frameAnchor: "top-left",
                  text: (d) => "Residual = " + d.meanResid
                })
      )
    ],
    facet: {
      data: summary_OB_OJS_all_T_selected_cohort_ind,
      x: "ind"
    }
  })

# Production model - river 4 (IS) Simple CJS model integrated with a growth in weight model to get phi, p, and growth estimates to develop a production model. In all the models below, 1 = not observed and 2 = observed in the input *encounter data*.\ Encounter data are available [here](https://drive.google.com/drive/folders/1UmPv49xL-mzUBndqjRISylNO3q3Zs-cm) in the `eh.csv` file. Weight data are in `weight.csv` ```{r} #| label: librariesModelsNimbleProduction #| echo: false library(getWBData) library(lubridate) library(kableExtra) library(GGally) library(nimble) library(nimbleEcology) library(MCMCvis) library(tidyverse) library(targets) library(shiny) library(tibble) ``` ```{r} #| label: globalModelsNimbleProduction #| include: false knitr::opts_chunk$set(warning = FALSE, message = FALSE) ``` Using model #3 from `modelsCMR_Growth_NN_OB.qmd` as a staring point for the models, but adapting the model by\ Prob not 1) Extending AgeInSamples from 1-11 to 1-x to allow bigger fish to be present for the production estimates. Prob not 2) Loop over first\[i\]:lastAIS so fish have the chance to survive to large size.\ 3) Add sample random effect structured by cohort to allow cohort effects on growth. ## By Cohort, modelNum 3: phi(i,t) \* g(i,t), p(i,t) model Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i) Probability of survival (phi) model structure: ``` logit(phi[i,t]) <- phiInt[i,t] + phiBeta[1,i,t] * weight[i,t] + phiBeta[2,i,t] * weight[i,t] * weight[i,t] + phiBetaCohort[cohort[i]] ``` Probability of capture (p) model structure: p(t,i) \<- pInt(t,i) Growth rate (gr) model structure: gr(t,i) \<- grInt(t) + cohort(i) Survival/growth rate interaction: Multiplicative Model code is in `./models/production/modelProduction_OB_functionsToSource.R`\ Model is run 'by hand' in `./models/modelProduction_OB_makeFile.R`\ Functions for this qmd file in `./models/qmdProduction_OB_functionsToSource.R`\ Model runs:\ ### How many ageInSamples to include? ```{r} #| label: production all <- tar_read(cdWB_electro_target) table(all %>% filter(river == "wb obear")|> dplyr::select(ageInSamples)) cohorts <- 2002:2014 ``` ### Retrieve model results ```{r} #| label: productionSource #| cache: false library(targets) # Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html # # # To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs source('./models/production/modelProduction_OB_functionsToSource.R') source('./models/production/qmdProduction_OB_functionsToSource.R') modelNum <- 3 # phi * growth ``` ### Model code ```{r} # all cohorts have the same code. show here load('./models/production/runsOut/production_OB_model_3_2002__mostRecent.RData') d$modelCode ``` ```{r} #| label: productionFunction getSummaryByCohort <- function(cohort, modelNum = 3) { # Load the data load(paste0('./models/production/runsOut/production_OB_model', '_', modelNum, '_', cohort, '__mostRecent.RData')) # Input data eh <- tar_read_raw(paste0('eh_OB_', cohort, '_target')) get_MCMCTrace_OB(d) MCMCplot(object = d$mcmc, params = c("phiIntTOut"), ref_ovl = TRUE) MCMCplot(object = d$mcmc, params = c("phiBetaT"), ref_ovl = TRUE) MCMCplot(object = d$mcmc, params = c("p"), ref_ovl = TRUE) MCMCplot(object = d$mcmc, params = c("grIntT"), ref_ovl = TRUE) # Create the summary list summary_list <- list( eh = eh, stats = kable(as_tibble(d$runData), caption = "Run statistics"), runTime = paste0('Run time = ', round(d$runTime, 3), ' ', attr(d$runTime, "units")), #traces = get_MCMCTrace_OB(d), summaryMod3_tt_growth = MCMCsummary(object = d$mcmc, params = c("phiIntTOut", "p", "grIntT", "phiBetaT"), round = 3) %>% rownames_to_column(var = "var"), summary_OB = get_summary_OB(d, eh) ) return(summary_list) } #s=summary_OB[, order(summary_OB$ind, summary_OB$t)] ``` ### Summaries by cohort ::: {.panel-tabset} #### 2002 ```{r} #| label: s2002 #| cache: TRUE s2002 = getSummaryByCohort(2002, modelNum) ``` #### 2003 ```{r} #| label: s2003 #| cache: TRUE s2003 = getSummaryByCohort(2003, modelNum) ``` #### 2004 ```{r} #| label: s2004 #| cache: TRUE s2004 = getSummaryByCohort(2004, modelNum) ``` #### 2005 ```{r} #| label: s2005 #| cache: TRUE s2005 = getSummaryByCohort(2005, modelNum) ``` #### 2006 ```{r} #| label: s2006 #| cache: TRUE s2006 = getSummaryByCohort(2006, modelNum) ``` #### 2007 ```{r} #| label: s2007 #| cache: TRUE s2007 = getSummaryByCohort(2007, modelNum) ``` #### 2008 ```{r} #| label: s2008 #| cache: TRUE s2008 = getSummaryByCohort(2008, modelNum) ``` #### 2009 ```{r} #| label: s2009 #| cache: TRUE s2009 = getSummaryByCohort(2009, modelNum) ``` #### 2010 ```{r} #| label: s2010 #| cache: TRUE s2010 = getSummaryByCohort(2010, modelNum) ``` #### 2011 ```{r} #| label: s2011 #| cache: TRUE s2011 = getSummaryByCohort(2011, modelNum) ``` #### 2012 ```{r} #| label: s2012 #| cache: TRUE s2012 = getSummaryByCohort(2012, modelNum) ``` #### 2013 ```{r} #| label: s2013 #| cache: TRUE s2013 = getSummaryByCohort(2013, modelNum) ``` #### 2014 ```{r} #| label: s2014 #| cache: TRUE s2014 = getSummaryByCohort(2014, modelNum) ``` ::: ### Combine cohort summaries ```{r} #| label: bindSummaries sAll <- tibble(.rows = 0) |> add_column(!!!s2002$summary_OB[0,]) # Fill tibble for (cohort in cohorts) { summary_obj <- get(paste0("s", cohort))$summary_OB summary_obj$cohort <- cohort sAll <- bind_rows(sAll, summary_obj) } sAll$cohort_ind <- paste0(sAll$cohort, "_", sAll$ind) ``` ```{r} #| label: ojsDefine #| cache: FALSE #ojs_define(numTags_OJS_mod3 = dim(s2002$eh$tags)[1]) # all cohorts have the same eh ojs_define(summary_OB_OJS_all = sAll) #ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3)) ``` ### Select cohort(s): ```{ojs} //| label: summary_OB_OJS_all_T summary_OB_OJS_all_T = transpose(summary_OB_OJS_all) ``` ```{ojs} //| label: selectCohort //| import { range } from "d3"; viewof selected_cohort = Inputs.select(d3.range(2002, 2015), { label: "Which cohort?", value: 2002, step: 1, multiple: true }) ``` ```{ojs} //| label: filterCohort summary_OB_OJS_all_T_selected_cohort = summary_OB_OJS_all_T.filter((d) => (selected_cohort.includes(d.cohort)) ) ``` ### Select one or more individuals: ```{ojs} //| label: selectInd uniqueInds = [...new Set(summary_OB_OJS_all_T_selected_cohort.map(d => d.cohort_ind))].sort(); viewof selectInd_mod3 = Inputs.select(uniqueInds, { label: "Which fish?", value: 1, step: 1, multiple: true }) summary_OB_OJS_all_T_selected_cohort_ind = summary_OB_OJS_all_T_selected_cohort.filter((d) => (selectInd_mod3.includes(d.cohort_ind)) ) ``` Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish's cohort. ### Plot survival Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish's cohort. ```{ojs} //| label: plotZ Plot.plot({ width: width, height: 350, inset: 10, color: { scheme: "greys" }, x: { label: "Age/season combination" }, y: { label: "Probability of survival" }, marks: [ Plot.frame(), Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind, { x: "t", y: "pSurv" }), Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind, { x: "t", y: "enc8", fill: "orange" }), Plot.line(summary_OB_OJS_all_T_selected_cohort_ind, { x: "t", y: "pSurv" }), Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind, { x: "first", y: 1, "stroke": "green" }), Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind, { x: "last", y: 1, "stroke": "red" }) , Plot.text(summary_OB_OJS_all_T_selected_cohort_ind, Plot.selectFirst({ x: 11, y: 1, text: d => d.cohort }) ) ], facet: { data: summary_OB_OJS_all_T_selected_cohort_ind, x: "ind" } }) ``` ```{ojs} //| label: plotWeight Plot.plot({ width: width, height: 350, inset: 10, color: { scheme: "greys" }, x: { label: "Age/season combination" }, y: { label: "Standardized body mass (mg)" }, marks: [ Plot.frame(), Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind, { x: "t", y: "mean_weight" }), Plot.line(summary_OB_OJS_all_T_selected_cohort_ind, { x: "t", y: "mean_weight" }),      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind, { x: "t", y: "weightDATA_std", fill: "orange" }), Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind, { x: "first", y: 3, "stroke": "green" }), Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind, { x: "last", y: 3, "stroke": "red" }), Plot.text(summary_OB_OJS_all_T_selected_cohort_ind, Plot.selectFirst({ x: 1, y: 4, frameAnchor: "top-left", text: (d) => "Cohort = " + d.cohort }) ), Plot.text(summary_OB_OJS_all_T_selected_cohort_ind, Plot.selectFirst({ x: 1, y: 3.5, frameAnchor: "top-left", text: (d) => "Residual = " + d.meanResid }) ) ], facet: { data: summary_OB_OJS_all_T_selected_cohort_ind, x: "ind" } }) ```